Resonance and propulsion performance of a heaving flexible wing 



Sebastien Michelin^' ^'j^] and Stefan G. Llewellyn Smith^ 

'Department of Mechanical and Aerospace Engineering, Jacobs School of Engineering, 
University of California San Diego, 9500 Gilman Drive, La Jolla CA 92093-0411. 
^Ecole Nationale Superieure des Mines de Paris, 
60-62 Boulevard Saint Michel, 15212. Paris Cedex 06, France. 
(Dated: June 15, 2009) 

The influence of the bending rigidity of a flexible heaving wing on its propulsive performance in 
a two-dimensional imposed parallel flow is investigated in the inviscid Hmit. Potential flow theory 
is used to describe the flow over the flapping wing. The vortical wake of the wing is accounted 
for by the shedding of point vortices with unsteady intensity from the wing's trailing edge. The 
trailing-edge flapping amplitude is shown to be maximal for a discrete set of values of the rigidity, 
at which a resonance occurs between the forcing frequency and a natural frequency of the system. 
A quantitative comparison of the position of these resonances with linear stability analysis results is 
presented. Such resonances induce maximum values of the mean developed thrust and power input. 
The flapping efflciency is also shown to be greatly enhanced by flexibility. 



I. INTRODUCTION 



Unlike terrestrial animals that can use solid friction and a fixed support, insects and fishes must generate from the 
surrounding fluid the lift and thrust forces necessary to their motion in their environment [TJ |2l O |4l |5] . Beyond the 
fundamental interest of understanding the mechanism of insect flight and fish swimming, recent research on propulsion 
in fluids has also been motivated by the development of micro-aviation vehicles (MAV) and of more efficient propulsion 
techniques based on biomimetics. 

Insects use thin flapping wings to generate an unsteady flow around them to produce these forces. The flow is 
characterized by a relatively large Reynolds number Re = UL/v with U the typical wing velocity, L its characteristic 
chord and v the kinematic viscosity of the surrounding fluid: Typically, Re ~ 100-5000 [5]; in that range, the forces 
on the wings are dominated by the pressure contribution and viscous effects are concentrated in thin boundary layers 
near the solid's boundary. These boundary layers separate during the unsteady wing motion and roll up into strong 
coherent vortices [B] that carry momentum away from the insect, thereby generating the propulsive forces. The flow 
around the insect is highly unsteady, and an ongoing research challenge resides in the ability to describe the forces on 
the flapping structure without solving explicitely for the details of the flow [TJ |H1 IHl HHJ IH] ■ 

An important physical insight into the generation of propulsive forces by flapping or deforming solids has been 
provided by the study of active motion and deformation. In such experimental [TH [13], theoretical [TJ [2] and 
numerical studies [TSJ [TBI [13 [H] , the position of the solid is prescribed and the influence of the swimming stroke on 
the propulsive performance is analyzed. 

The recent development of experimental imaging techniques has however shown that most insect wings are not 
purely rigid and that their deformation is not entirely controlled by the animal: the wing can experience large passive 
deformations during the stroke period [19^ under the action of the outside flow and its internal bending rigidy, whose 
spanwise and chordwise distribution results from the venation pattern of the wing |20l I21j . The forcing is generally 
applied by the insect on its wing through a main axis whose rigidity is signiflcantly higher than the rest of the 
structure (e.g. the leading edge of the wing). An important challenge now is to understand the impact of such passive 
deformation on the propulsive performance of a flapping structure. In particular, one may be interested in potential 
reduction of energy usage induced by the flexibility, and examine if the values of the rigidity for these structures lie 
within an "optimal" range for which the flapping efficiency is the highest. The definition of optimality in this work 
in terms only of thrust production and propulsive efficiency is purposely restrictive: from a biological point of view 
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many other factors must be taken into account to define the optimal structure of an insect wing, including but not 
limited to material resistance and manoeuverability. 

The purpose of the present study is to investigate numerically the effect of flexibility on the propulsive characteristics 
of a flapping appendage. In the following, this structure will be referred to as a wing, understanding that the model 
could also be applied to a flsh flu if it is allowed to deform passively under the effect of the flow and of its bending 
rigidity. Solving for the coupled motion of a flexible solid and a fluid is computationally challenging and expensive, 
primarily due to the coupling occurring on a moving boundary whose position is a priori unknown and must be solved 
for as well. Popular techniques to overcome this difficulty are the use of coupled fluid and solid solvers using fltted 
grids [22] and immersed boundary methods |23i |24J . The use of low-order models for the flexible wing also simplifies 
the computation while still retaining important physical results on the reaction of the body to the fluid flow |25l [26] . 

The model used here focuses on a simplified two-dimensional propulsion problem using a fiexible wing of infinite 
span, actuated at its leading edge in a purely heaving motion and reacting passively to the fiow forces and its internal 
elasticity. The present model does not aim to represent a particular fiying or swimming pattern, but rather considers 
a one-degree-of- freedom forcing to focus on the influence of flexibility on the performance of the apparatus. In the 
limit of high Re, viscous forces are neglected and the viscosity's influence on the flow is retained in this potential flow 
formulation by the irreversible shedding of vorticity from the trailing edge of the flapping structure, in the form of 
point vortices whose unsteady intensity is determined so as to satisfy the regularity condition at the solid's trailing 
edge [23128]. 

A similar approach was recently proposed by Alben [29^ for a pitching elastic sheet using a vortex sheet represen- 
tation of the wake. Optimal values of the solid's rigidity were discussed in the limit of negligible solid inertia (that is 
particularly relevant in the case of fish swimming) and of infinitesimally small displacements of the solid. The present 
work builds on these results and considers the general case of non-linear deformations of the sheet with non-negligible 
inertia (as is the case for a flapping insect wing). The use of the unsteady point vortex model rather than the full 
vortex sheet description also allows for a simpler treatment. Similar resonance patterns are observed and a theoretical 
argument is provided for their origin and position. The relation between thrust or drag production and vortex wake 
structure is also investigated. 

In section [IT] the fluid-solid model is presented and the propulsive performance quantities of interest are defined. 
Section III discusses briefly the numerical methods used as well as the existence of a periodic regime. Section IV then 
investigates the influence of the solid's rigidity on the propulsion forces and efficiency and relates them to the structure 
of the solid's wake. Peaks of thrust are observed for particular values of the rigidity and in section |V| the occurrence 
of such peaks is showed to correspond to a resonance between the forcing frequency and the natural frequencies of the 
fluid-solid system. Finally, section VI presents some general conclusions and discusses the limitations of the model. 



II. DESCRIPTION OF THE MODEL 



A. Solid model 

The following two-dimensional model for the flapping structure is considered (see Fig. [T| . The wing is represented 
by an elastic sheet of chord L and inflnite span, clamped at its leading edge on an attachment pole of negligible 
thickness, actuated by the operator (e.g. main body of the insect). The operator applies a purely vertical motion to 
the sheet's leading edge, whose orientation is constrained to be stricly horizontal. The vertical position and orientation 
of the wing at the leading edge are then: 

h{t) = /io(l - cosujt) and 6'o(t) = 0, (1) 

so that A — 2ho and / — lu/2tt are respectively the amplitude and frequency of the flapping motion (to avoid confusion, 
uj or its non-dimensional form will be referred to as the angular frequency in the following). The flapping wing has 
a chordwise flexural rigidity per unit length B and a mass per unit area ps- Its thickness is negligible compared 
to L. The wing is placed in a uniform horizontal flow J7oo of density p. The motion of the wing's leading edge is 
entirely prescribed by (jlj but the rest of the wing has a purely passive motion in response to its internal elasticity, 
the leading-edge forcing and the pressure forces applied by the surrounding flow. 

In the following, L, Uoc and p are used as reference quantities to non-dimensionalize the problem. The properties 
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FIG. 1: Heaving flexible wing in a steady axial flow. The heaving motion of amplitude A is imposed at the leading edge and 
vortices are shed from the trailing edge. 



of the elastic sheet are characterized by the mass ratio /i and non-dimensional rigidity 77 defined as 

Ps B 

and the leading-edge forcing is characterized by the non-dimensional forcing amplitude e and frequency / 

'-"i-T. -'-it- 

The forcing Strouhal number is defined in accordance with previous experimental studies |12j as 

St={^^2ef. (4) 

The motion of the wing is described using an inextensible Euler-BernouUi beam representation [2S] . We are interested 
in large displacements of the wing so all non-linear geometric terms must be included. The linear Euler-BernouUi 
assumption remains valid if the curvature radius of the wing is much larger than its thickness, which is assumed here. 
The position of the wing is described using complex notation as ^(s, t) — x(s, t) + iy(s, t) and its orientation is defined 
as 9{s,t) with < s < 1 the curvilinear coordinate along the wing. The classical notation for the complex flow 
velocity is also used here: w = u — iv, with {u,v) the cartesian components of the velocity vector. The conservation 
of momentum for each element of the wing and the inextensibility condition can be written as 

mC= [(T-ir;0..)e'«]^-i[p]±e'^ C. = e'^ (5) 

where [p]^ is the pressure difference between the top and bottom sides of the wing and T is the wing tension that 
must be solved for at each time-step to enforce the inextensibility condition. The clamped-free boundary conditions 
imposed by the forcing (jlj are 

C(0,i) = Co(t) = ie[l-cos(27r/i)], 0(O,t) = O, (6) 

e,{i,t)^esAi,t)^T{i,t) = o. (7) 

Equations (|5])-(|7]) can be rewritten as a system for 9 and T only [2H]: 

T,, - e^^T = -[pfe, - 2rj0,9,,s - vOl ~ (8) 
^ie = - riOssss + {T+ riODOss + 2Tse, (9) 

e{Q,t) = e,{\,t) = ess{i,t) = T{\,t) = o (lo) 

Ko + M^ / e'''(i^-^2jds'ds=-r(O)+i770,,(O)-i^ [pj^e'^ds. (11) 
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B. Representation of the flow around the flapping wing 



The flow around the wing is taken potential. The expected boundary layer separation at the trailing edge and its 
subsequent roll-up into vortices is represented in this inviscid formulation by a discrete shedding of point vortices 
from the trailing edge Ce — Ci^i^) of the flapping wing (Fig. [T]). The intensity of the last shed vortex is determined so 
as to cancel exactly the square-root singularity arising in the velocity field due to the presence of the flat sharp corner 
[22l EHl EHl [31] . When the intensity of the unsteady vortex reaches a maximum, a new vortex is created from the 
generating corner, thereby expressing the irreversible nature of vortex roll-up in the formalism of this inviscid model. 
From that point on, the intensity of the previous vortex is frozen. The unsteady point vortices satisfy the modified 
equation of motion 



Zn + {Zn - Ce)^^ = Wn 
^ n 



(12) 



where z„ and r„ respectively refer to the position and intensity of the point vortex. Wn is the desingularized complex 



velocity at the vortex position and the overbar denotes a complex conjugate. Equation (12) is known as the Brown 



Michael equation [32 , and enforces the conservation of fluid momentum around the vortex and associated branch 
cut in an integral sense [27\. The omission of the corrective term on the left-hand-side would lead to an unphysical 
unbalanced force on the branch cut linking the vortex to its generating corner |33j . 

The shedding of vorticity from the leading edge is neglected here, as we focus mostly on situations where the angle 
of attack remains small at the leading edge [29] . Alternatively, this representation can also be seen as the limit case 
of a smoothed leading edge of very small curvature radius (as in an airfoil profile for example). Only one unsteady 
vortex is shed at a time (from the trailing edge). Noting N the number of vortices at a particular time, all vortex 
intensities r„ are therefore independent of time except for the last one Tj^^t). 

In the absence of viscosity, the tangential velocity of the flow can be discontinuous across the wing. The potential 
flow around the wing is then computed by representing the infinitely thin solid as a bound-vorticity distribution k 
[m [Ml [35] . The complex fiow velocity is obtained by superposition of the flow at infinity and the contribution of the 
bound and wake vorticity 



w(z, t) 



1 
27ri 



Kds 



N 

E 



(13) 



The regularity condition at the trailing edge imposes w{(f.,t) 7^ oo or equivalently K{l,t) = 0. The total circulation 
at infinity is conserved and, assuming the system is started from rest, must be zero at all time. The bound-vorticity 
K is the solution of a singular Fredholm equation obtained by applying the continuity of normal velocity on the wing 
[18j . The corresponding system of equations for k and Fjv is 



1 

2^ 



Re 



.Wis) 



C(so) - C{s) 



K{s)ds = Im 



1 
2^1 



N 
n=l ^ 



-c 



„1 N 

/ K(s)ds + Vr 

"'0 , 1 



K{l,t) 



(14) 

(15) 
(16) 



From K, the pressure jump [p]^ across the wing can be computed by integration of Bernoulli's theorem along the 



wmg: 



/•So 

[p]='=(so) = / k{s)ds + k{so)wp{sq) 
Jo 

with Wp the principal value of the relative tangential velocity on the wing [28 

1 f K{s)ds ^ iF 



Wp{so) =Re 



,ie(so) 



ZTTiJo C(so)-C(s) 



-C(=^o) 



(17) 



(18) 
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C. Energy conservation 



From dsl , the conservation of energy can be written for the wing 



where 



Ek 



1 



\C\'ds, 



Ep) = Wp+V,n, 



(C(0)e-"'( 



ds, 



(r(o) + i770,,(o)) -r^me,{Q). 



(19) 

(20) 

(21) 

(22) 
(23) 



are respectively the kinetic and elastic potential energy of the wing, the rate of work of the pressure forces on the 
wing and the rate of work of the force and torque applied by the attachment pole on the rest of the wing. In the 
particular case of a purely heaving motion considered here, 9{Q) — and ^(0) — ih is purely imaginary, so 



V^n=^lOss{Q)h. 



(24) 



D. Propulsive performance 

We are interested in the thrust generated by the flapping wing. The forces applied on the leading-edge attachment 
are: 

— the elastic forces applied by the sheet on its attachment [T(0) — i770ss(O)] e'^'-^-', 

— the force applied by the operator (or animal) to prescribe the leading-edge motion F^p = F^^ 



the suction force created at the leading edge by the inverse-square root behavior of the pressure. This suction 
force is the limit of the suction force obtained on a smoothed contour when the curvature radius of the airfoil's 
leading edge tends to zero. This suction force is equal to ^29| |33 157] 



,ie(o) 



(lini Vs(l - s)k(s) 



(25) 



Neglecting the inertia of the attachment and defining the thrust (counted positively to the left) as T = —F^^, the 
force balance along the horizontal direction together with ([6])-([7]) leads to 

2 



r 



^(lim[y^(T^«:(s)]) -r(0) 



The instantaneous power input V by the operator is 

= Re (FopCo) = Fy^hit) = 770,,(O, t)h{t) = V^-. 



(26) 



(27) 



and is equal to the rate of work Vin of the attachment pole on the wing. Note that this equality would not hold if 
the motion of the leading edge were a combination of both heaving and pitching, as the suction force Fg would have 
a non-zero rate of work along the vertical direction. 

The useful power output is simply the rate of work TUco of the thrust force in the horizontal motion. In non- 
dimensional units, the fiapping efficiency is then defined as the ratio of the average developed thrust to the average 
input power 



(28) 
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where (.) is the averaging operator over a flapping period 

{g) = - r 9it)dt, (29) 

T Jo 

with T = 1// the non-dimensional period of the flapping motion and the positive part of V. In the following, 
the mean power input is understood as (V^), thereby assuming that the animal cannot store and reuse the energy 
possibly extracted from the fluid if < during a fraction of the flapping period. For the range of parameter values 
used here, it was observed that V > for most of the flapping period and (V^) ^ (V), except for very rigid wings. 
The results and discussions presented here are therefore not affected by this choice. 

Finally, the following (non-dimensionalized) quantities are deflned for convenience: 

— the trailing-edge peak-to-peak flapping amplitude V, 

— the intensity of the wake F^, defined as the mean value of the amplitude of the successive vortices (positive and 
negative) , 

— the induced velocity of the wake vortices V defined as the horizontal velocity of the wake vortices relative to 
the imposed unit flow. V is positive if the wake vortices move faster than the background flow, and negative 
otherwise. 



III. NUMERICAL SIMULATION OF THE INITIAL VALUE PROBLEM AND CONVERGENCE TO A 

PERIODIC STEADY STATE 

Equations (|8|-(12l and (14)-(17l are solved numerically expanding 6{s,t) into a flnite series of Chebyshev polyno- 
mials of the flrst kind and using a semi-implicit second-order time-stepping scheme. Taking advantage of the linear 
relation between [p]^ and k, added inertia terms can be isolated from the part of the pressure that can be explicitely 
computed at each time step, thereby avoiding the use of an iterative solver and greatly enhancing the computational 
efliciency \28^ . 

The system is started from rest. At t = 0, the horizontal flow is ramped up to its long time unit value and the 
motion of the leading edge ([ij is imposed. After a transient regime of a few heaving periods, a permanent periodic 
regime is achieved for large enough values of the rigidity rj (Fig. ^f). 

However, the harmonic heaving forcing can lead to highly unsteady behaviors if rj becomes too small: below a 
certain critical value of the solid's rigidity rjmil^), the purely passive elastic sheet (as in the flag problem) becomes 
unstable to fluttering modes and flapping can occur even in the absence of leading-edge forcing |28} 1351 138j (f?m(/i) 
was found equal to rjm — 2 10"^ for fj, — 0.2 and ri,n — 4.8 10~^ for fj, = 2, using the same point vortex model |1H])- In 
such cases, the spectrum of the trailing-edge motion can display several peaks, corresponding to the forcing frequency 
and to the frequency of the unstable modes (Fig. |2|d). The motion is of large enough amplitude for the regime to be 
non-linear and mode coupling is also expected. As rj is reduced further, the periodicity is lost and the power spectrum 
is full; in such a case, determining averaged quantities is not possible anymore. 

This explains the difficulty to observe a steady permanent regime when rj is decreased below the critical value r/m- 
This difficulty was not present in the linear study by Alben \^ as the inertia of the solid was neglected (^ ^ 0) . The 
inertia of the solid is essential to the development of fluttering instability and in the limit ^ 0, all the modes are 
linearly stable [251 1351 |3H] . In the following, unless indicated otherwise, the range of parameter values is chosen such 
that a periodic state is achieved. 



IV. WING FLEXIBILITY AND PROPULSIVE PERFORMANCE 



In this section, the behavior of the propulsive performance (mean thrust, mean power input and efliciency) is 
studied when the rigidity rj of the wing is varied. Several values of the forcing frequency /, forcing amplitude s and 
mass ratio fj, were investigated. 
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(a) /i = 2 and rj = 0.2 
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(b) ^ = 2 and 77 = 0.04 



FIG. 2: (Top) Frequency spectrum and (Bottom) time evolution of the trailing-edge absolute vertical displacement for fJ, = 2, 
e = 0.05 and / — 5/2n. (Left) rj = 0.2 > r|m{^J■ = 2) lies in the stability region for the purely passive elastic sheet. The 
power spectrum displays only one peak (thick arrow) at an angular frequency of u/ = 2nf = 5. (Right) The rigidity rj is 
below its critical value r;m(pi = 2) = 0.048 and the elastic sheet is unstable to fluttering. The power spectrum displays two 
main peaks: one corresponding to the forcing frequency (w/ =5, thick arrow) and one corresponding to the unstable fluttering 
mode {iUr ~ 2.3, thin arrow) that matches the flapping frequency observed in the purely passive case e = 0. In both cases, 
small peaks can be seen for oj ~ 3a;/ = 15, corresponding to the third harmonic of the forced flapping. The fact that only odd 
harmonics appear in the tail motion was already previously observed in the case of a passive flag [28) . 



A. Optimal flexibility for thrust generation and propulsion efficiency 

For given mass ratio, forcing amplitude and frequency, the mean thrust, mean power input and propulsive efficiency 
were computed for each value of the rigidity 77. As a general result, starting from the case of a rigid wing (77 — s- 00), 
the mean thrust (T) and power input (7^"^) both increase when flexibility is introduced in the problem; however, the 
former increases faster and the resulting efficiency is an increasing function of flexibility (decreasing function of rj) 
for large values of 77 (Fig. [s]). When the solid's flexibility is further reduced, (T) and {V'^) display successive peaks 
occuring for the same values of rj. The propulsive efficiency r displays in general one large peak (in general for a 
different value of the rigidity rj than the thrust peaks) before dropping sharply to zero as the mean thrust vanishes 
(drag-thrust transition). The increase in the flapping efficiency and developed mean thrust for a flexible wing is 
significant compared to the case of a rigid wing: for / = 5/27r, the peak value of the mean thrust can be greater than 
twice its value in the rigid case and the efficiency can increase from 27% to almost 60% (Fig. ^jp). Similar behavior 
is observed for / = I/tt (Fig. [s^). Flexibility has therefore a significant impact on the performance of the propulsive 
apparatus considered here. 

Figure [3] shows the evolution of the propulsive properties for the lowest possible value of rj leading to a permanent 
periodic regime. In the higher frequency case (/ = 5/27r), the drag-thrust transition was clearly observed while at 
lower frequency (/ — I/tt), this transition is less well defined as unsteady phenomena start developing around the 
same value of the rigidity. Several factors lead to unsteadiness of the problem, including transitions in the wake 



behind the flapping wing and development of the fluttering instability for low values of the rigidity (see section III I 



The behavior of (T), (7^^) and r collapse rather well for leading-edge flapping amplitudes up to 50% of the wing's 
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FIG. 3: Evolution of the mean thrust (T) (top), power input {V^) (center) and propulsive efficiency (bottom) with the wing's 
rigidity t), for a flapping wing of mass ratio /i = 0.2, flapping frequency (a) / = l/vr and (b) / = 5/27r and forcing amplitude 
e = 0.1 (solid-circle), e — 0.2 (dashed-star) and e — 0.5 (dotted-square). For comparison purposes, the mean power input and 
thrust have been normalized using their rigid case value as 77 — > oo. 



length and lower frequencies (Fig. |3p). However, for higher frequencies, the influence of the forcing amplitude can be 
seen for values of e as low as 0.1 (Fig. [3|d). As a general result, the propulsive efficiency and the normalized mean 
thrust and power input are observed to decrease with e when all other parameters are held fixed. The optimal values 
of 77 for maximal thrust or maximal efficiency are also observed to increase with e. 

It must be emphasized here that the decrease of the achievable thrust and power input with the forcing amplitude 
e are only relative to the rigid case: for higher forcing amplitudes, the absolute mean thrust and power input are 
larger in magnitude and, in the rigid limit, are observed to scale like e^/^. 



B. Wake structure and thrust production 

1. Evolution of the wake structure with increasing flexibility 

For given mass ratio, frequency and forcing amplitude £, the variation of the flexibility of the wing induces important 
changes in the developed thrust and energy use. As rj is decreased from the rigid case {rj 00), the wake behind 
the flapping wing also undergoes important modifications. Figure |4] shows the evolution of the fiow pattern around 
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(a) Rigid case rj = 890 
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FIG. 4: Streamlines of the flow over the flapping wing for fj. — 0.2, e = 0.1 and / — 5/2tt and decreasing rigidity ij. The positive 
(resp. negative) vortices are represented with upward- (resp. downward-) pointing triangles with sizes scaled to the intensity 
of the vortices. The streamlines are plotted for t = 20 when the leading edge crosses the horizontal axis, and the leading-edge 
forcing for < t < 20 was the same in all three cases. 



the heaving wing for varying rj. In the high-rigidity case (Fig. |4p,), the trailing-edge deflection is small. The wake 
vortices are arranged in a reverse Von Karman vortex street, in which vortices with positive (resp. negative) intensity 
are positioned above (resp. below) the horizontal axis. This arrangement induces an acceleration of the fluid on the 
horizontal axis to the right and the vortices have a higher velocity than the imposed background flow (positive V). 
This jet carries momentum to the right which is consistent with the formation of a mean thrust on the wing. The 
formation of a reversed Von Karman vortex street in thrust-producing propulsion schemes is well known and has been 
observed in several experimental studies [12 [T5] . 

When rj is decreased, the solid is more flexible and the trailing-edge flapping amplitude V increases. As a result, 
the intensity of the wake vortices increases with decreasing rigidity and so does the thrust generated by the flapping 
motion. Figure Qb) shows the streamlines and arrangement of the wake vortices for rj = 3.2 which corresponds to the 
flrst peak in thrust production on Fig.jsjb). The phase between the leading and trailing-edge displacements has also 
been modifled compared to the rigid case. The increase in the vortex advection velocity can be seen as the distance 
between two successive vortices is slightly increased (the vortex shedding frequency is unchanged and equal to the 
forcing frequency). 

As a comparison. Fig. ^c) shows the flow pattern for the value of 77 corresponding to the thrust-drag transition. 
One observes immediately that the intensity of the wake vortices has decreased significantly and the width of the 
vortex street has almost vanished. In intermediate Re experiments, the formation of a classical Von-Karman street 
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FIG. 5: (Square) Induced wake vortices velocity V behind the flapping wing for varying -q and ^ = 0.2, e — 0.1 and / — 5/2tt. 
The predicted value obtained for an infinite reversed Von-Karman vortex street [39] is plotted for comparison (solid line). 



is associated with drag on the generating sohd body. It appears here that, even at the transition between thrust 
and drag production, the reversed Von Karman pattern persists, although largely weakened. This result is consistent 
with the observation in recent experiments |13j that the thrust-drag transition does not necessarily occur at the same 
time as the transition in the wake structure. A theoretical argument for this difference is also presented in section 



IV B 2 Furthermore, in purely inviscid simulations, the Von Karman street in the drag-producing case is generally 



much weaker and instead the vortices tend to align along the axis 



2. Classical and reversed Von Karman streets and thrust/drag production 

Here, theoretical results on vortex streets are used to understand the relationship between drag/thrust production 
and vortex wake structures. We focus on highly periodic cases where the structure of the reversed Von-Karman or 
classical Von Karman streets is easily identified. 

For a given vortex wake intensity and vortex arrangement, the advection velocity of the vortices is the superposition 
of the incoming flow velocity (equal to 1 in non-dimensional units) and of the induced velocity V of the vortex street, 
which is itself a direct function of the width of the vortex street (vertical distance between the two rows of opposite 
sign vortices), its intensity and its wavelength (horizontal distance between two successive vortices of identical sign). 
If the vortex street were infinite, its induced velocity V would be |39j 

V^^S^tanhf^V (30) 



2a \ a 



where r^, is the magnitude of the vortices, b the width of the vortex street and a its horizontal wavelength. In (30 1, 
the following convention is chosen: b > (resp. b < 0) for a reversed (resp. classical) Von-Karman street in which 
positive vortices are located above (resp. below) the horizontal axis. If the shedding frequency (equal to the fiapping 



frequency /) is fixed, then a — (1 + V)//. Measuring the width of the vortex street b and its intensity F^u, (31 1 leads 
to a non- linear equation for Vl, the predicted value of the induced velocity, that can be solved numerically for given 
b and F^u 

Figure |5] shows very good agreement between the measured value V and the expected value Vl for varying rj (the 
other parameters taking the same values as in Fig. |4|, particularly above the thrust-drag transition at 77 = 8 10~^ 
(Fig. |5|. This agreement shows that the induced velocity is mostly determined by the neighboring vortices, and the 
semi-infinite or infinite nature of the vortex street does not infiuence the induced velocity significantly. 

The drag on a solid body placed in a uniform parallel fiow and shedding vortices in the pattern of a staggered vortex 
street (regular or reversed Von-Karman street) was computed by Von Karman using the conservation of momentum 
around the solid body and part of its wake |36[ I40j . The only important elements in the derivation are again the 
intensity T^, wavelength a and width b of the vortex street. In particular, the motion (other than the main translation) 
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FIG. 6: (Square) Mean thrust (T) produced by the flapping wing for varying i] and /i = 0.2, e — 0.1 and / — 5/2n. The 
predicted value 7s (32 1 based on impulse conservation [361 140j is also plotted for comparison (solid line). Ts was computed 



from ( 32 1 using the values of the vortex street intensity, width and induced velocity obtained with our model. 



and deformation of the body are irrelevant in Von Karman's derivation. These results can easily be generalized to 
the case of a reversed Von-Karman street and in our notation, the predicted mean thrust Ts is obtained as 

T5 = -;P^ + ^(1 + 2V), with a = (1 + V)//. (32) 

This theoretical prediction compares very well to the results obtained for the mean thrust in our simulations (Fig.|6|. 
In particular, the transition between thrust and drag production is well reproduced. The agreement is lost at low 
values of rj: For such low values of the rigidity, the highly regular structure of the wake is also lost as some natural 
modes of the passive elastic sheet become unstable to fluttering and the motion of the wing loses its strong periodicity, 
making the definition of a, and b difficult. 

In ( p2| ), the thrust consists of two terms. The second term is positive in the case of a reversed Von Karman street 
(6 > 0) and negative in the classical Von Karman street (with 6 < and V < 0, as long as V > —1/2, which always 
occur in the weak Von Karman streets observed here) . The presence of the first term (which is always negative and 
therefore always leads to drag production) is responsible for the experimentally- and numerically-observed difference 
between the thrust/drag transition and the transition in the wake structure. As in the passive fiapping fiag case 



[281 135] , it is possible to have net drag produced by a weak reversed Von Karman street as the first term in ( 32 1 
dominates the second one. 



C. Resonance and optimal flexibility for thrust production 

In Fig.|3|b), successive peaks in the mean thrust created by the heaving wing can be observed, but also a peak in 
the mean drag (or negative thrust) below the drag-thrust transition that occurs around 77 ~ 8 . 10~^. Figure [t] shows 
the evolution of the mean thrust, fiapping amplitude V, wake intensity and wake induced velocity V for the same 
values of the parameters /i, e and / as in Fig. [3j 

A very clear correlation is observed between the occurence of the maxima (in magnitude) for the mean thrust (or 
drag) and for the other quantities, which confirms the argument presented in the previous section: the mean thrust 
and drag peaks are created by an increase in the fiapping motion at the trailing edge, where the vortex wake is 
formed. An increase in V (with a constant fiapping frequency) induces higher relative velocity at the trailing edge, 
and therefore stronger shed vortices. While this increase in r„ with V is physical, it is however not possible to find 
a simple scaling of r„i with V, as other factors that depend on 77 must be taken into account (e.g. the orientation of 
the trailing edge relative to its velocity). The induced velocity on the vortex street V is therefore also increased and 
the wake carries a larger fiuid momentum downstream, thereby creating a greater thrust on the profile. 

It is also interesting to notice that the maximum drag observed for 77 ^ 5. 10^^ is also associated with a maximum 
in D and r,„. In that case, the large amplitude of motion at the trailing edge opposes the imposed fiow and creates 
a net drag on the body. 



12 




FIG. 7: Evolution of the mean thrust (T) (solid), trailing-edge flapping amplitude O (dashed-circle) , vortex wake intensity Fm 
(dotted-square) and induced vortex velocity V (dotted-triangle) for fi = 0.2, e = 0.1 and / = 5/2n. All quantities have been 
normalized by their rigid-case value (r; — + oo). 



V. INFLUENCE OF FLEXIBILITY ON THE FLAPPING AMPLITUDE AND MODE SHAPE 

A. Resonances between the forcing frequency and the natural frequencies of the system 

In this section, we are interested in the origin of the maxima in the traihng-edge flapping amphtude V. The 
successive peaks in I? as 77 is varied suggest a resonance phenomenon. By varying the rigidity of the sohd 77, the 
natural frequencies of the system are also modified. For an elastic sheet in vacuum, these frequencies scale like 
^ B I ps- Looking for solutions of (|8])-( 11 ) in the linear limit with no fluid forcing, the fundamental angular frequencies 



LOQn — 27r/„ of a clamped-free elastic sheet in vacuum are obtained in non-dimensional form as 

ujQn = -^ni/ with 1 + cosh A„ cos A„ = 0. (33) 



The ratio \fiijr\ can be thought of as the non-dimensional time-scale associated with the frequency of the sheet's 
natural oscillations in vacuum, or alternatively as the outside flow velocity non-dimensionalized by the characteristic 
velocity associated with the sheet's properties. 



1. Influence of fi and e on the position of the resonance peaks in T) 



However, the resonance peaks observed in Fig. [7] do not correspond to a resonance with the natural frequency 
of the elastic sheet in vacuum. If this were the case, then the resonance would be achieved for all values of the 
mass ratio n at the same value of \/rj/ p. Such a coincidence does not occur (Fig. [s]): the position of the successive 
resonances is actually strongly influenced by the fluid-solid inertia ratio ji. This difference makes sense physically, 
as the eigenfrequencies of the system are modified by the presence of the forcing horizontal flow, and such effects as 
added inertia are expected to be important. Instead, the natural frequencies of the system {wing -t- outside uniform 



fiow} should be considered. A comparison with linear analysis prediction is proposed below in section VA3 



Before comparing the numerical results to the linear analysis, we study the influence of the forcing amplitude e 
on the position of the resonances. It was noticed in section [IV A| that the normalized thrust and power input follow 
similar patterns for values of e up to 0.2 to 0.5 depending on the value of p. In an attempt to determine the exact 
value of rj leading to a resonance, we study how the resonance peaks in T) are modified as e is varied between 0.01 and 
0.5. Starting from small £, one observes that the increase in e induces a shift of the resonance peaks toward larger 
values of rj and a smoothing of the resonance peaks (Fig. [9|. Also, the convergence toward the limit case of small 
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FIG. 8: Evolution of the normalized trailing-edge flapping amplitude V with yjrj/ ^ for e = 0.05, / = 5/27r and for ^ = 0.2 
(solid-squares) and /j, = 2 (dashed-circle). The flapping amplitude was normalized using the asymptotic rigid case limit {rj oo). 




FIG. 9: Influence of the forcing amplitude e on the position of the resonances in the trailing-edge flapping amplitude for /i = 0.2 
and / — 5/2n. The results are plotted for e = 0.01 (solid-circles), e = 0.05 (dashed), e = 0.1 (dash-dotted), e = 0.2 (dotted) 
and e = 0.5 (solid). The trailing-edge flapping amplitude was normalized by its rigid-case value {-q oo). 



£ is faster for the second and third resonance peaks than for the first one (labeling the peaks from the right as 77 is 
decreased from the rigid case limit). 



2. Absolute and relative trailing-edge flapping amplitude 



In the previous sections, T) was defined as the peak-to-peak amplitude of the trailing edge fiapping motion in the 
laboratory frame (thereafter referred to as absolute flapping amplitude). In the following, we are also interested in the 
motion of the wing in the frame moving with the leading edge. In this frame, the peak-to-peak trailing edge flapping 
amplitude is V* (thereafter referred to as relative flapping amplitude). V and V* are not necessarily equal because 



of the phase difference between the motion of the leading and trailing edges (Fig. 10). This delay is, in general, a 
decreasing function of the solid's rigidity. In the limit of a rigid wing rj — > 00, the trailing edge flaps in phase with the 
leading edge with the same amplitude {V* = and T) — 2e). As rj is decreased, a delay appears between the motion of 
the leading and trailing edges: the elasticity of the wing takes more time to carry along its length the signal imposed 
at the leading edge. As a result, the relative amplitude V* increases. 

One of the main consequences of the appearance of such a delay is the non-coincidence of the first resonance peak 
(with largest rj) in the fiapping amplitude, whether T> or V* is considered. The position of the subsequent peaks are 
not significantly affected (Fig. llOk) 
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FIG. 10: (a) Comparison between the absolute flapping amplitude T> (solid) and relative flapping amplitude T>* (dashed) for 
the trailing edge, for /i = 0.2, / = 5/27r and e = 0.1. (b) Phase difi'erence between the leading-edge and trailing-edge motions. 



3. Comparison with the natural frequencies of the system 

In this section, the conjecture that the traihng-edge flapping amphtude peaks are due to a resonance with the 
natural frequencies of the {wing + imposed flow} system is tested by comparing the value of the rigidity 77 leading 
to a peak value for T) to the value of rj for which the forcing frequency matches a natural frequency of the system. 
The position of the peaks in P is measured in the limit of small forcing amplitude (typically e = 0.01). The natural 
frequencies of a purely passive wing clamped at its leading edge in a uniform flow are determined using the linear 
stability analysis method developed by Kornecki |41j . A brief summary of this method is given in the Appendix [] 
For given /z and 77, the eigenfrequencies of the system are computed. The lowest frequencies are also associated with 
the modes of lowest order (those with the longest wavelength). Here, the following equivalent problem is considered: 
for a given /i and /, we want to find the values 77 for which a mode of the passive elastic sheet in axial flow has the 
particular frequency /, regardless of its growth rate. 



Figure 11 shows the position in the (77, /)-plane of the first resonances observed in the trailing-edge flapping ampli- 
tude T) for small heaving amplitude e, starting from the rigid case 77 — > 00. Although e is small, T) can be significant 
(greather than 20e at the resonance in the case of the heavier wing (/i = 2), thereby representing more than 20% of 
the wing's length). These results are compared to the natural frequencies of the system {wing + parallel fiow} as 
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FIG. 11: Position of the resonances for the trailing-edge flapping amphtude in the (77, /)-plane obtained using the present 
model (symbols) with e = 0.01 and (a) /i = 0.2 and (b) \i = 2. The different symbols correspond to the nature of the mode: 
modes 1 have no neck in the flapping enveloppe (circles); modes 2 (squares), modes 3 (upward-pointing triangle) and modes 4 
(downward-pointing triangles) have respectively 1, 2 and 3 necks in their motion envelope. The black symbols correspond to 
resonances in the absolute flapping amplitude T). Open symbols correspond to resonances in the relative flapping amplitude D* 
when they differ from the resonances in T). The position of the resonances is compared to the prediction of the linear analysis 
(dashed) for the natural frequency of the purely passive elastic sheet (or flag) in axial flow (see Appendix [l. 



predicted by the linear analysis and very good agreement is found between the numerical results and the theoretical 
predictions. The position of the resonances in T)* are also indicated. Resonances in V* and T) coincide except for the 
first resonance peak in the case of the smaller mass ratio \i. 

The linear analysis seems to underpredict slightly the values of 77 corresponding to the resonance for a given 
frequency This difference is consistent with the amplitude of discrepancy between the point vortex model and the 
linear stability analysis observed in the study of a purely passive elastic sheet or flag [2S] . Two other factors can also 
explain the small discrepancy in the results: 

— As pointed out, the motion of the wing is not infinitesimally small here, even for e — 0.01; it is therefore 
possible that non-linear effects modify the exact position of the resonances. In the previous section, the effect 
of increasing e was shown to shift the resonance peaks toward larger 77, particularly for the first peak, and this 
could account for a significant part of the observed discrepancy. 
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The amplitude of flapping in the laboratory frame is considered here. However, even for small e, this amplitude 



differs from the relative flapping amplitude defined in the frame moving with the leading edge (see Fig. 10 1, 
because of the existence of a non-zero phase between the motion of the leading and trailing edges. It was 
observed in Fig. [lO]that the resonance in relative amplitude occurs for smaller values of ij, particularly for the 
first resonance (largest rj). In the linear analysis, both amplitudes (relative and absolute) are identical since the 
leading edge is held fixed. For ^ = 0.2, the agreement is improved for the position of the first resonance if V* 



is considered instead of V (Fig. 11 1 



Figure 11 also indicates the nature of the observed mode, in particular the number of necks in the envelope of the 
wing's motion. One observes that resonances located on a same branch of the linear analysis prediction share the 
same general structure, and the number of necks is consistent with that predicted by the linear analysis: for given fi 
and 77, the lowest frequency mode has the longest wavelength and no neck in its envelope. The next lowest frequency 
corresponds to a mode with one neck, and so on for the successive frequencies. The evolution of the mode shape for 
varying 77 is studied in more detail in section |VB[ 



B. Evolution of the flapping mode shape with the flexibility of the profile 

To confirm that the maxima of flapping amplitude actually correspond to resonances between the forcing frequency 
and the natural frequency of the passive elastic sheet in a parallel flow, the evolution of the flapping mode shape with 
r] is now considered. For comparison with the case of a purely passive flexible elastic sheet, the mode shape is defined 



as the envelope of the motion of the wing in the frame moving with the leading edge. Figure 12 shows the mode shape 



in permanent regime at the values of 77 leading to peak values of the absolute trailing-edge flapping amplitude and 



to the minima between two successive peaks. Note that on Fig. 12 case B seems to correspond to a wider envelope 
than case A, although the absolute flapping amplitude is smaller for B than for A. This is the result of the change 
of frame: if one considers the flapping amplitude in the moving frame, the position of the maxima differs from peaks 



in absolute flapping amplitude (see Fig. 10 1. Cases A, C and E correspond to the resonances while B, D and F 
correspond to local minima of the absolute flapping amplitude. 

The mode shapes C and E are structurally similar to the envelope of the first two flapping modes observed for the 
passive flexible flag which is consistent with C and E corresponding to resonances between the flapping frequency 
and the flapping modes 2 and 3 (one- and two- neck modes respectively). A corresponds to mode 1 of the passively 
flapping elastic sheet, which was not observed in the case of the flapping flag study [55] as it is always stable [32 and 
therefore does not lead to spontaneous large-scale flapping. 



VI. CONCLUSIONS 



Using a reduced-order model for the flow past a two-dimensional heaving flexible wing, the influence of the wing's 
flexibility on its propulsive performance (mean thrust, flapping efficiency) was investigated. Starting from the purely 
rigid case, we observed that the flexibility of the wing allowed for a larger trailing-edge flapping amplitude, thereby 
generating a stronger wake and an increased mean thrust. The energy usage also increases with the introduction of 
flexibility, but more slowly than the mean thrust, resulting in a net increase in the flapping efficiency with reduced 
rigidity. This efficiency gain can be significant (up to twice the efficiency of the rigid case) . While the mean thrust and 
power input display several peaks when the rigidity rj is decreased from the purely rigid case, the flapping efficiency 
displays one wide peak before falling sharply as 77 nears the value leading to the thrust-drag transition. Below this 
threshold, the wing is too flexible to communicate momentum to the flow and instead starts creating a net drag on 
the leading-edge attachment. 

The relationship between thrust production and wake structure was then investigated, taking advantage of the 
discrete representation of the wake. Analytical predictions for the induced vortex street velocity and mean thrust in 
terms of the vortex street strength and spatial arrangement were successfully compared to our simulation results. 

The peaks in mean thrust were found to correspond to maximum values in the trailing-edge amplitude, and shown 
to be the result of the resonance between the forcing frequency of the heaving motion and the natural frequencies 
of the system. A quantitative comparison showed very good agreement between the optimal values of the solid's 
rigidity and the linear analysis predictions for the resonances position. The existence of these resonance phenomena 
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FIG. 12: Evolution of the mode shape with rj for fj, = 0.2, e = 0.1 and / — 5/2n. (Top) Traihng-edge flapping amplitude in 
the stationary frame T> (solid) and in the frame attached to the leading edge V (dotted). Note that the value of -q leading to 
a resonance in T> does not necessarily correspond to the value of rj leading to a resonance in V* . (Bottom) Mode shape plotted 
for the value of rj indicated on the top panel. The position of the wing in the frame moving with the leadinge edge is plotted 
every = 0.06. 



was further confirmed by comparing the flapping mode shape to the mode shape observed for a freely flapping elastic 
sheet (e.g. the flag problem [55]). 

The natural frequencies of the system are strongly dependent on both the solid's flexibility and the ratio of fluid 
and solid inertia, and so are the optimal values of the solid's rigidity. For a slender neutrally buoyant fish fin, the 
inertia ratio /i can generally be neglected [55]. However, in the case of an insect wing, the small thickness-to-chord 
ratio is balanced by the large difference in density for the fiuid and solid, and the present analysis shows that the 
mass ratio fi plays an important role in defining the optimal value of the flexibility. 

The model presented here used a potential flow representation and the shedding of point vortices to describe the 
highly-unsteady flow around the insect wing. Vortices were shed only from the trailing edge. This approximation 
is reasonable for small angles of attack (proportional to the Strouhal number St in the case of a purely heaving 
motion). Then, the vorticity shed at the leading edge is negligible or merges with the vorticity shed by separation of 
the boundary layers at the trailing edge [T5], leading to the shedding of two individual vortices every flapping period 
and so-called 2S wakes [43] . However, when the flapping amplitude and frequency are increased, more complex wakes 
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are expected as vorticity is shed from both the traihng and leading edges, inducing for example the formation of the 
so-called 2P wake, where a vortex pair is shed during each half period. The representation of the leading-edge vortex 
falls however beyond the scope of inviscid methods such as the present point vortex model or a vortex sheet approach 
as the effect of viscosity can not be neglected [23144]: the leading-edge vortex is expected to remain close to the body 
for a sufficiently long time to interact with the boundary layers. The use of point vortices also restricts this method to 
two-dimensional problems and therefore does not allow the study of such effects as wing-tip vortices that are expected 
to influence significantly the flight performance. The present study however does not aim at reproducing the exact 
flow around an insect wing but to present a test case where the effect of flexibility can be isolated from other factors 
such as leading-edge and wing-tip vortices. 

Despite these limitations, the present method offers the advantage of a considerable reduction in computational 
cost for two-dimensional fluid-solid simulations. This is particularly attractive for situations where the cost of full 
numerical simulation is prohibitive and a large number of simulations are required (typically in the case of optimization 
problems). 

The present work was purposely limited to a one-degree of freedom flapping pattern (pure heaving) in order to limit 
the number of free parameters and focus on the fundamental effect of solid flexibility on the flapping performance. 
In future work, more realistic flapping schemes should be considered to follow more closely the flapping pattern of 
an insect wing. In particular, a combination of heaving and pitching should be considered to determine the influence 
of the relative phase between heaving and pitching motions on the results presented in this work. The interaction 
between two flapping flexible sheets should also be investigated to understand the lift and thrust generation by insects 
with multiple pairs of flexible wings (e.g. dragonfly j45]) or the efficiency of fish schooling, and to complement recent 
experimental studies on multiple passive fiexible filaments [ISl ITZl H51 
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APPENDIX: NATURAL FREQUENCIES OF A PASSIVE FLEXIBLE SHEET IN AXIAL FLOW 

We present here a brief summary of the linear stability analysis of a clamped-free fiexible sheet or fiag using a vortex 
sheet representation of the wake as developed by Kornecki jH]. More details on the method and the full calculation 
can be found in Refs. [HJ [SD]. Small vertical displacements of the sheet h{s,t) are considered (0 < s < 1) with 
I ft. I <C 1 so that the equation of motion of the solid ^ becomes in linearized form 

fih = -rihssss - ^P- (A.l) 

The wake is represented by a continuous distribution of vorticity 7(5, i) along the horizontal axis (s > 1) advected by 
the fiow (this continuous shedding thereby differs from the discrete approach presented in the present paper, but is 
better suited to linear stability analysis). The pressure forcing can be decomposed into two parts jUJISniEI]: a non- 
circulatory part due to the flow created by the solid's motion with no net circulation around the solid and a circulatory 
part due to the flow created by the vortex wake (the absence of circulation at inflnity requiring the existence of a net 
circulation around the solid equal to the opposite of the wake vorticity). Considering a decomposition of the solid's 
motion onto normal modes of the form h{s,t) — Re [y(s)e"^'], y{s,t) satisfles j41j 



19 



s 



2iw 

TT 



log 



1 



7rv/a;(l - a;) Jo x-^ 
1 ' 



1-C 



(icjy + y^)d^ 



{iujy + y') 
2a:- 1 



- x) V X 



1 - X 



C{u) 



where C{ijj) is the Theodorssen function [50] : 
r(2) 



2ij(^) (I) 



(2) 



(A.2) 



(A.3) 



(A.4) 



and Hij '{x) = Jy(x) — iY^{x) {v = 0, 1) are Hankel functions of the second kind fS^^. 

The eigenvalue problem (A.2 1 for w is solved numerically using a Galerkin method: y{s) is decomposed along the 
first TV eigenmodes of the clamped-free beam in vacuum y{s) = ^ a„^("^(s) with ^^"^^s) satisfying 

V.(")(o) = ^(")(o) = = ^^(1) - 0, 

and A„ are the successive positive roots of 1 + cos A„ cosh A„ — 0. ( |A.2 ) is replaced by the nonlinear eigenvalue 
problem 
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FIG. 13: For fi — 0.2, evolution of (Left) the growth rate and (Center) frequency of the first four modes (of lowest positive 
frequency): mode 1 (solid), mode 2 (dashed), mode 3 (dash-dotted) and mode 4 (dotted). (Right) Critical stability curve as 
obtained in the (/i,77)-plane from the linear stability analysis. 



For a given N, these matrices can be precomputed. Then, for given /i and rj, the eigenvalue problem (A. 5) is solved 
iteratively using a Newton-Kantorovitch algorithm 
frequency w/27r for fi - 
below r]jn = 2. 10" 



Figure 13 shows the evolution of the growth rate -~lm{uj) and 



0.2 and varying rigidity rj. In the rigid case (large rj), all modes are stable. As rj is decreased 
one or more modes become unstable (mode 3 then mode 4). The critical stability curve separates 
the region of the {fi,rj)-pla,ne where the elastic sheet's state of rest is stable and the region of the parameter space 
where at least one mode is unstable (Fig. 13). 
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